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P5 Abstract 



(/3 A method that uses fuzzy clustering algorithms to achieve particle identification based on pulse shape analysis is presented. The 
O fuzzy c-means clustering algorithm is used to compute mean (principal) pulse shapes induced by different particle species in an 
c/3 automatic and unsupervised fashion from a mixed set of data. A discrimination amplitude is proposed using these principal pulse 
^^ shapes to identify the originating particle species of a detector pulse. Since this method does not make any assumptions about 



Oh' 



the specific features of the pulse shapes, it is very generic and suitable for multiple types of detectors. The method is applied 
to discriminate between photon- and proton-induced signals in CsI(Tl) scintillator detectors and the results are compared to the 
well-known integration method. 
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1. Introduction 



■<!::;j- ' In many experimental setups, particle identification (PID) 

^^ is an important task. A well-known approach to identify the 

^^ species of a particle hitting a detector is the use of pulse shape 

. . analysis (PSA), which exploits that the detector response de- 

K>" pends on the particle species. 

k>( ■ A very widely-used method for achieving particle identifi- 
'V^ cation using PSA is the integration method (see e.g. UJ-lsl]), 
Ci which uses the ratio of the integrals over two regions of the 
detector signal to distinguish between different particle species. 
A property of this method is that the regions that provide the 
best distinction have to be determined beforehand and — in or- 
der to achieve the best results — separately for each individual 
detector As detectors age and are exposed to radiation, their 
responses might change over time. Changes of external pa- 
rameters (e.g. temperature) might also affect the pulse shape 
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properties and require a recalibration of parameters. Thus, this 
calibration procedure, which requires human supervision, has 
to be repeated regularly. 

For very large setups with thousands of detectors, e.g. the 
CALIFA barrel ||1 of the R^B experiment at FAIPQ, this is 
hardly feasible. This is also true for detectors used in environ- 
mental field inspection with strongly varying conditions. For 
this reason, it is highly desirable to develop a procedure that re- 
quires the least human supervision possible to perform the PID 
calibration of a detector system. 

We have recently shown that fuzzy clustering algorithms can 
be used to perform the calibration for y-neutron discrimination 
in a fully automatic and unsupervised manner for liquid scintil- 
lator detectors ||5|]. Due to the success of this method, we apply 
the same prescription in this work to the y-proton discrimina- 
tion using scintillating CsI(Tl) crystals, which have been shown 
to have distinct scintillation characteristics for different incident 
particles JQ^. We furthermore derive a means of using the 
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information obtained from the clustering to perform the PID 
which is suitable for use in a signal processing environment, 
e.g. for online-identification using fast digital electronics. In 
addition, the proposed procedure is independent of the partic- 
ular features of the pulse shapes and thus represents a general 
approach suitable not only for CsI(Tl) detectors, but also for 
other detector types. 

This paper is structured as follows: section ll.il gives a short 
introduction on the topic of clustering algorithms and outlines 
the algorithm on which this manuscript concentrates. In sec- 
tionlj, the experimental setup is described during which the set 
of signals was recorded that we use to demonstrate the method. 
Footnote |5| focuses on the preprocessing of the raw signals that 
is employed in order to get stable and correct clustering results. 
The results obtained from the clustering process are then used in 
the prescription described in section |4l There, we also summa- 
rize the well-known integration method and present a different 
approach, based on a specific distance measure, that does not 
need the results from the clustering process. 

1.1. The fuzzy c-means algorithm 

Clustering algorithms provide a means of assigning a set of 
A^ vectors in L-dimensional space to C groups, called clusters. 
This assignment is done based on a notion of the distance be- 
tween a vector from the set and the centroid of a cluster The 
cluster centroids themselves are not known beforehand, but are 
determined by the algorithm. 

The degree of membership of the /-th vector to the j-\h clus- 
ter is designated by a membership value Ujj, which is normal- 
ized so that the total degree of membership ^^^j m,j = 1 for 
every vector Based on the range of the membership values, 
clustering algorithms may be distinguished into two categories: 

• boolean clustering algorithms with m,j € {0, 1}, i.e. a vec- 
tor belongs exactly to one cluster. 



degree of membership of any vector should be equal to unity. 
Here, Xj denotes the /-th vector of the set, Cj is the centroid of 
the j-th cluster and || ■ || is the Euclidean norm; A^ and C denote 
the number of vectors and clusters, respectively. Note that C is 
a parameter of the algorithm that has to be known beforehand. 

The parameter m e [1, oo) describes the fuzziness of the al- 
gorithm: For m - \, the algorithm degenerates to the boolean 
k-means clustering algorithm, and for the limit m -^ oo, all 
membership values become m^ - I/C, yielding a completely 
fuzzy clustering. In the range of m = 1.1 to m - 2.0 the ob- 
served dependence on m of the clustering results is very weak 
and we use m = 1 .2, which yielded the best results, in the fol- 
lowing discussion. For a more detailed discussion of the m- 
dependence we refer to section 6.1 of [|5|]. 

Employing Lagrange multipliers, one can derive the neces- 
sary conditions for a minimum of 7,„, viz. 
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and find a local minimum by iterating these equations starting 
from randomly chosen M,y obeying the constraints in ([T]i. The 
random choice of the m^ ensures an unbiased application of the 
clustering procedure. 

In the following, we extensively use the isomorphism be- 
tween a sampled detector signal s{k), where I < k < L de- 
notes the sample index, and a vector s represented as an L- 
tuple s = (^(fc))/(:eii,L] with respect to an arbitrarily chosen or- 
thonormal Cartesian basis of K^. We nevertheless keep the 
notation s(k) to emphasize that these signals are functions of 
(discretized) time. 



fuzzy clustering algorithms that use the full range of uij e 2. Experimental setup 



[0, 1] and thus can assign a single vector to multiple clus- 
ters with different degrees of membership. 

In the following we will focus on fuzzy clustering algorithms 
because they comprise the more general case. 

The concrete fuzzy clustering algorithm we successfully em- 
ployed earlier |5] is the fuzzy c-means algorithm |)9nlOJ. It can 
be posed as a minimization problem 

A' C 



min Jm^YjJ] "O'll^' ~ ^ill^ 
subject to 

N 

2m,;>0, ;■£{!,. ..,C} 

c 
^My = 1, ie{l,...,N} , 



(1) 



The data used to demonstrate the method proposed in this 
work was taken during a beam time at the Mayer-Leibnitz lab- 
oratory (MLL) in Garching near Munich. The MLL tandem 
proton accelerator (see e.g. OlI - llSll ') provided a proton beam 
with energy of 24 MeV. The protons impinged on a 70 (xm CD2 
foil target. The elastically and inelastically scattered protons, as 
well as photons stemming from nuclear reactions inside the foil, 
were detected by a 8x4 array of CsI(Tl) detectors iH. [l4l Hsl 
at an angle of 37° with respect to the beam axis. The aver- 
age count rate of a single crystal was of the order of 200 Hz 
resulting in a very small amount of pile-up signals (see also be- 
low). The crystals were read out using Hamamatsu S8664-1010 
APDo A Mesytec MPRB-16 shaping pre-amplifier, which 
also provided the bias voltage for the APDs, was used to obtain 
the detector signal. This signal was sampled by a 60 MHz 12-bit 
analog-to-digital converter on a FEBEX2 FPGAO board. The 



i.e. minimizing the total intra-cluster distance Jm subject to the 
constraints that no cluster should be empty and that the total 
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Multi Branch System data acquisition software Ill6i] recorded 
per-event signal traces, each 2000 samples in length, and addi- 
tional data generated by the FPGA board for offline analysis. 

3. Signal preparation 

To ensure a stable clustering process, the signal set is pre- 
pared by performing baseline subtraction, noise reduction, time 
alignment and removal of distorted signals. These steps are de- 
scribed in the following, the procedure is the same as in the 
application of the method to liquid scintillators ISl] . 

First, the signal baseline voltage is estimated by taking an av- 
erage over the first 100 samples of the trace, where no signal is 
present, and then subtracted from the signal. Next, the signals 
are filtered twice using a FIRo low -pass filter with a cutofl" fre- 
quency of 1.8 MHz. This is done to reduce noise superimposed 
on the signal which could impede the clustering. The cutoff 
frequency of the low -pass filter is large with respect to the char- 
acteristic time scale of the signal: the shortest rise times are 
larger than 1 ^s. Thus, distortions introduced by this step are 
minimal. 

For the clustering to work, the signals have to be aligned in 
time. This alignment is performed using a digital CFE|f|: First, 
a search for the signal maximum r^nax = ''(^max) is performed, 
with A:max denoting the sample number where the maximum is 
assumed. Starting at k^^ax and going to earlier times (smaller k), 
the index ^^fd of the first sample with value less than 0.2 r^^x is 
determined. In case of signal pileup, the location of the global 
signal maximum ^max is most likely the maximum of the sec- 
ond (piled-up) pulse and the search for the constant fraction 
point yields a sample number that is far away from the signal 
maximum. The constant fraction point, however, will in most 
cases reside on the rising edge of the first pulse. Thus, the signal 
maximum is determined again in the region [A;cfd-50, A:cfd+450] 
and a new constant fraction search is performed. If the dififer- 
ence between the old and new constant fraction points is larger 
than 50 samples, the new point is taken as fccfd and the process 
is repeated. That way, only the first of the overlapping pulses is 
selected. Note that a sub-sample alignment as performed in yl] 
is not necessary because the sampling period is small compared 
to the rise times of the signals. 

The energy of the event is then estimated by taking the aver- 
age of the five samples at k^ax - 2, ... , k^-^x + 2. Because the 
clustering results deteriorate when the ratio between photon- 
and proton-induced signals is too far from unity ||5|, the energy 
range between 2MeVee and lOMeVge, which contains a com- 
parable number of signals from both species, is selected (see 
fig. 111). All signals whose energies lie outside this region are re- 
moved from the set. Since only the proton pulse shape shows an 
energy dependence, it would be possible to extend the cluster- 
ing to higher energies by first applying the fuzzy c-means to the 
selected region which contains both types of signals and then 
running the algorithm in the region with E > lOMeVee only 




Fig. 1: Energy deposited in tlie detector plotted against the distance Z)(.v, ■) 
to the mean signal shape of the full set, given by the D distance measure (see 
section l4!3l for details). The measure separates both particle species well. In the 
region above 8 MeVee only proton signals are present, whereas photon signals 
dominate in the region below 4MeVee. The (shaded) region between 2MeVee 
and 10 MeVee was selected because it contains a sufficient number of signals of 
each signal class. 



for the proton cluster while keeping the photon cluster centroid 
fixed. 

The signals are then cropped so that only 50 samples before 
and 450 samples after k^fd are retained. This ensures that the 
constant fraction point is the 5P' sample in all cropped signals 
Si(k). The cropped signals are normalized such that the signal 
energy is equal to unity. 

Because a contamination of the signal set with distorted sig- 
nals (e.g. due to pile-up or clipping) may lead to inferior clus- 
tering results, these signals have to be eliminated (see [5] for 
details). Therefore, a mean signal shape 



1 ^ 

m = - ^ si(k) 



(4) 



finite impulse response 
'constant fraction discriminator 



is calculated, where A^ is the number of signals still in the set. 
All signals with a euclidean distance larger than 2 .7 to this mean 
shape are considered distorted and eliminated from the set. This 
treatment was sufficient for the signal set at hand, removing 
only about 0.5% of the signals. For sets with a larger number 
of distorted signals a more sophisticated method, such as hier- 
archical clustering, may be necessary. 

After performing the correction of the CFD timings, pileup 
did not play a significant role for the present signal set. Any 
signal that contains a piled-up pulse starting inside the cropped 
region is considered distorted and discarded by the last step. 
Piled-up pulses starting after the cropped region are ignored. 

4. Particle identification 

In the following sections, two methods for achieving particle 
identification are presented, the first of which is the well-known 
integration method. The second method uses fuzzy clustering 
to determine a principal pulse shape for each particle species 
and may be combined with various discrimination functions 
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Fig. 2: (color online) Typical normalized single signal shapes for protons and 
photons of a Csl(Tl) detector with APD readout passed through a shaping am- 
plifier The shaded areas mark the integration Umits used by the integration 
method. 
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Fig. 4: Energy deposited in the detector plotted against the Q[/Q2 ratio used 
by the integration method to discriminate particle species. The proton signal 
shapes show an energy dependence in the Q1/Q2 ratio while the photon shapes 
are largely independent of energy. 
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the particle species that induced it. The ratio of two integrals 
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Fig. 3: (color online) Q1/Q2 histograms for two different detectors. The peaks 
are well-separated, but their positions in both histograms differ slightly due to 
dift'erent characteristics of the pulse shapes for the two detectors. 



that make use of the principal pulse shape. Two discrimina- 
tion functions are presented and their performance is compared 
to the performance of the integration method. 

The comparison is done using the signals of two neighboring 
detector crystals near the center of the 8x4 array. We show 
a PID histogram for each of the investigated methods and each 
detector and, in order to provide a benchmark, we give a Figure- 
of-Merit (FOM) in the respective figures. The FOM is defined 
as 



M 



O"! +0-2 



(5) 



where yUi,/Z2 and 1x1,0-2 are the means and standard deviations 
of gaussian functions fitted to the two peaks of each histogram. 
The FOM is thus a mesure of how well the respective method 
separates the two particle species for the data set at hand. 

4.1. Integration method 

Because dift'erent particle species induce distinct processes in 
the scintillator, the signal shapes will diff'er and thus the integral 
over a part of the signal will give different values depending on 



61 ^ ULa, ^W 

22 zL,m 



(6) 



over two dift'erent regions [ai,bi] and [02, b2] is taken to make 
the value independent of the total integral and thus usable for 
particle identification. A histogram of the Q1/Q2 ratios for the 
set of signals is then created. Each particle species in the signal 
set produces a peak in the histogram at a characteristic Q1/Q2 
ratio and cuts can be applied to perform particle identification. 
In general, the integration limits have to be adjusted to the char- 
acteristics of the detector signals to get good identification re- 
sults. 

The integration limits for the signal set at hand are illustrated 
in fig. |2t The Q1JQ2 histogram calculated using these limits is 
depicted in fig. |3| for two different detectors. The peaks in the 
two histograms have different positions for the indiviual detec- 
tors (due to slightly different characteristics of the pulse shapes) 
and these positions cannot be known beforehand. They have to 
be determined by using a fitting procedure or by manual inter- 
vention. The FOM values for the two detectors are also given in 
fig.|3l Due to the wide proton peak in the histogram for detector 
A the value is significantly smaller 

Figure |4J shows the signal energy as a function of the Q1/Q2 
ratio for the signals in the present set. The photon signal shape 
is largely independent of energy. However, the shape of the 
proton signals shows an energy dependence, causing an energy- 
dependent change in the Q1/Q2 ratio. 

4.2. Fuzzy clustering 

The fuzzy clustering method achieves particle identification 
by comparing detector signals to prototype signals for each par- 
ticle species. These prototype signals do not have to be con- 
structed manually, but can be derived from a mixed set con- 
taining a large number of signals originating from the different 
species. This makes the fuzzy clustering method a two-step 
procedure, with the derivation of the prototype signals as first 



■^ 

¥ 




-0.4 



50 100 150 200 250 300 350 400 450 500 

k 

Fig. 5: Weighting functions obtained from l|8) for two detectors, used for the 
calculation of the amplitude histograms depicted in fig. [3 
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Fig. 6: (color online) Histograms of discrimination values for two detectors. 
The histograms are calculated from the same signal sets as the Q1IQ2 ratios 
from fig. Q The proton and photon peaks are situated around a discrimination 
value of and 1 respectively and the proton and photon signals form completely 
disjunct distiibutions. 



and the assignment of individual signals to particle species as 
second step. 

The first step is performed by applying the fuzzy c-means 
clustering algorithm to the signal set. The clustering yields 
a prototype signal Cj{k), given by the determined cluster cen- 
troids, for each particle species and a matrix of membership 
values Uij denoting the degree of membership in each cluster 
for each signal. The intra-cluster variance per sample k 



1 '^ 



(7) 



with the total cluster size Uj - 2, m,j, and other statistical quan- 
tities can easily be calculated from these results. 

The membership values can be used themselves to perform 
particle identification. However, they tend to values very close 
to unity or zero when the fuzziness parameter m is small and 
are thus hardly comparable to the results of other particle iden- 
tification methods. 

It is thus desirable to derive a different means for achieving 
particle identification from the results of the clustering process. 
The information obtained from the clustering method facilitate 
the use of an improved integration method. Instead of setting 
the integration regions manually, the prototype signals can be 
used to determine a set of weighting functions 



Wijik) = 



Ci(k) - Cj(k) 
A]{k) + Ajik) 



(8) 



These weights are largest where the prototype signals c,(A:) of 
the clusters under consideration differ the most with respect to 
their variances Af{k). As there is considerable scatter inside 
each cluster, the widths of both signal distributions have to be 
taken into account. This is achieved by relating the difference 
between the prototype signals to the sum of the intra-cluster 
variances of both clusters. The weights obtained for the present 
data set are shown in fig.lJ. 

Using the weights Wij(k), a discrimination value for each pair 
of particle species 



A,j[s(k)]^Y,'',j(k)s(k). 

k=0 



(9) 



with L being the length of the signal in samples, can be assigned 
to each signal from the set. This step can be expressed as a 
scalar product of two vectors and is suitable for implementation 
on fast digital electronics. 

The discrimination value A/j for the two prototype signals 
can also be calculated. This is advantageous, because the dis- 
crimination values can be normalized so that a value of zero and 
unity is assigned to both prototype signals respectively: 



Aij[s(k)] = 



AijMk)]-Aij[cj(k)] 
Aij[ci(k)]-Aij[cj{k)] 



(10) 



With the help of this information, the peak positions in the dis- 
crimination value histogram corresponding to the two clusters 



are A,-^ - and Ay = 1 by definition. The widths of the peaks 
can be estimated using the intra-cluster variance of the discrim- 
ination values 

1 ^ 
AA2.(0 = 77 X •^ni(Au[s„(k)] - Aij[c,(k)]f , (11) 



u, 



which describes the variance of the peak produced by cluster 
/ in the histogram of discrimination values computed with re- 
spect to clusters /, j. These initial guesses can be used to per- 
form an automated least-squares fit to the histogram, which 
should converge rapidly due to the accuracy of the initial 
guesses. Then, cuts can be automatically derived from the fitted 
distributions. 

Figure |6| shows the histograms of discrimination values for 
the same signal sets as in fig. [j The peaks corresponding to 
protons and photons are situated near zero and one, as expected. 
The peak positions are also identical for both detectors, even 
though the pulse shapes for the two detectors differ slightly, 
leading to different distributions in the Q1/Q2 histogram (see 
fig. |3|). Compared to the integration method also the FOM is 
slightly improved (compare the given values in the respective 
figures). This shows that our approach of defining the weighting 
function without prior knowledge of the properties of the signal 
shapes gives reasonable results. 

The prescription for obtaining the discrimination amplitudes 
can be interpreted as a special case of Fisher's linear discrim- 
inant for vanishing covariance between the signal samples de- 
scribed in [ 10]. Although that assumption is not strictly true due 
to the initial low -pass filtering of the signals, which induces sig- 
nificant correlations between neighboring samples, the results 
provide an excellent means for discriminating between the par- 
ticle species in an unsupervised way. 

A similar idea has been pursued earlier for photon-alpha dis- 
crimination in [8]: their method D consists of a weighted in- 
tegration over the signal to obtain a value that can be used for 
discrimination. The difference to the procedure proposed in this 
work is how the weights are obtained: by averaging over two 
separate sets of known photon and alpha signals, the authors 
obtain two prototype signals. The difference between these two 
signals is then used to construct the weight by fitting a polyno- 
mial to the rising edge of the difference and attaching an expo- 
nential tail to it. The tail decay constant is then optimized for 
maximum cluster separation. 

Using fuzzy clustering to obtain the prototype signals renders 
the need for a separate calibration set for each particle species 
unnecessary. Fisher's linear discriminant ensures optimal dis- 
crimination independent of the specific form of the difference 
signal. 

4.3. Distance approach 

Instead of performing a weighted integration over the signal, 
a notion of the distance between two signals can be used to per- 
form particle identification. Interpreting the signals as vectors 
in a high-dimensional vector space, the most straight-forward 
way of defining a metric is using the euclidean distance between 
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Fig. 7: (color online) The scaled distance D to the proton prototype signal for 
the same data as in fig. Q and fig. @. Using the value of D for the photon 
prototype signal the distance can be scaled to result in and 1 respectively. 



the two signal vectors 



d(su S2) = ^J(sl - S2V(si - S2) . 



(12) 



This naive approach, however, has some drawbacks, because it 
does not take the properties of real signals into account. The 
signals contain noise which is accumulated when calculating 
the euclidean distance and may significantly increase the dis- 
tance between two otherwise very similar signals, producing 
false results. 

A way to circumvent this has been employed in jsl]: a dis- 
tance measure 



£'(*l,*2) = J](«l(^)-.V2(fc)) 



(13) 



is defined, which allows for noise contributions to cancel be- 
tween different samples. This distance measure can then be 
used to calculate a discrimination amplitude based on the dis- 
tance to one of the cluster prototype signals. Again, the ampli- 
tudes of the other prototype signals can be calculated and used 
to scale the discrimination amplitude so that the prototypes are 
at zero and unity, respectively. Figure |7| shows a histogram of 
the distances to the proton cluster centroid D(cp, ■), scaled to ob- 
tain a plot similar to the discrimination value histogram shown 
in fig.y. Compared to the method using the weighting function 
the determined values of the FOM for the distance approach are 
slightly worse and comparable to the integration method. This 
is probably not surprising since the weighting function explic- 
itly concentrates on the region of maximum distance between 
the clusters, taking statistical fluctuations into account. 

We want to mention in this context one pecuUarity of the D 
measure, which is that the distances between three points are 
tightly related: 

DisuS2) = J](sm-S2(k)) 

k 
= Yjisiik) - S3(k) + S3(k) - S2(k)) 



^D(suS3) + D{s3,S2) 



(14) 



A change of the reference signal for the calculation of the dis- 
crimination amplitude thus results only in a (irrelevant) constant 
shift of the amplitudes. For the two-species case this can be 
exploited by taking the overall mean signal (HJi as a reference, 
effectively removing the need for performing the clustering. In 
addition to that, no care has to be taken that the ratio of the sig- 
nals induced by each of the particle species under consideration 
is not too far from unity. The drawback is that the peak posi- 
tions and widths are not known beforehand anymore. We used 
this approach in fig. [ij in order to select the energy region over 
which the clustering is performed. 

5. Conclusions and outlook 

We propose a robust and reliable method to perform PID cal- 
ibration of a detector in an unsupervised manner It is applied 
to a data set of mixed photon- and proton-induced signals taken 
with CsI(Tl) detectors at moderate count rates and similar pop- 
ulation of proton and photon signals. This method uses fuzzy 
clustering algorithms to obtain prototype pulse shapes for each 
particle species in a mixed set of signals, which are used to de- 
rive a discrimination amplitude for each signal. The calculation 
of the discrimination values is not computationally expensive, 
i.e. once the weighting functions have been determined for a 
given detector they can be used to calculate the discrimination 
amplitudes for further data online using fast digital electronics. 
By recording detector signals for a small percentage of events 
and regularly performing clustering on this set, a change of the 
prototype pulse shapes (e.g. caused by a change of conditions) 
can be detected and acted upon by updating the discrimination 
amplitude. 

The amplitudes derived this way are used to assign signals 
to particle species using discrimination regions that are — in 
contrast to other methods — detector-independent. Furthermore, 
the results show that using the derived weighting functions the 
Figure-of-Merit, a measure of the separation between particle 
species, is comparable to the integration method. 

In an earlier work ||5[ we successfully applied the method 
to y-neutron discrimination in liquid scintillator detectors. In 
this work, we show that the same prescription can be used 
to achieve particle identification for y-proton discrimination in 
CsI(Tl) scintillators. This shows that the proposed method is 
indeed a generic one, which can be applied to many different 
detector types. The method is in principle also applicable for 
more than two particle species. This case should be investi- 
gated in upcoming works. 
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